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Abstract. We consider a simple Markovian class of the stochastic Wilson-Cowan 
type models of neuronal network dynamics, which incorporates stochastic delay caused 
by the existence of a refractory period of neurons. From the point of view of the 
dynamics of the individual elements, we are dealing with a network of non-Markovian 
stochastic two-state oscillators with memory which are coupled globally in a mean-field 
fashion. This interrelation of a higher-dimensional Markovian and lower-dimensional 
non-Markovian dynamics is discussed in its relevance to the general problem of the 
network dynamics of complex elements possessing memory. The simplest model of 
this class is provided by a three-state Markovian neuron with one refractory state, 
which causes firing delay with an exponentially decaying memory within the two-state 
reduced model. This basic model is used to study critical avalanche dynamics (the noise 
sustained criticality) in a balanced feedforward network consisting of the excitatory and 
inhibitory neurons. Such avalanches emerge due to the network size dependent noise 
(mesoscopic noise). Numerical simulations reveal an intermediate power law in the 
distribution of avalanche sizes with the critical exponent around —1.16. We show that 
this power law is robust upon a variation of the refractory time over several orders 
of magnitude. However, the avalanche time distribution is biexponential. It does not 
reflect any genuine power law dependence. 
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1. Introduction 

Network complexity pervades biology and medicine [I], and the human organism can be 
considered as an integrated complex network of different physiological systems [2j such 
as circulatory and respiratory systems, visual system, digestive and endocrine systems, 
etc., which are coordinated by autonomic and central nervous systems including the 
brain. The dynamics of the sleep-wake transitions during the sleep of humans and other 
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mammals presents one of important examples of such complex dynamics. In turn, 
functioning of the human brain presents is essence a network activity of the coupled and 
interrelated neurons surrounded by glia cells. The immense complexity of this subject 
matter [6] does not exclude, but rather invites thinking in terms of simple physical 
modeling approaches, see e.g. in Ref. [3], since even very simple physical models can 
display very complex behavior. The models of critical dynamical phenomena such as 
self-organized criticality (SOC) P0IHJ are especially important in this respect [3115] . 
Physical modeling can help to discriminate the physical and biological complexity from 
the complexity of mental processes, the “form within” [6], which is mediated but 
not determined in fine features by the background physical processes. The recently 
discovered complexity of the critical brain dynamics mm in essence does not have 
anything in common with the complexity of mental processes as it is already displayed by 
organotypic networks of neurons formed by cortical slices on a multi-electrode array [9]. 
Such physical complexity is in essence the complexity of crude matter that got self- 
organized following the physical laws. It thus belongs to statistical physics or system 
biophysics. Physical models such as SOC are especially important and helpful here. 

The Wilson and Cowan model ca presents one of the well-established models 
of neuronal network dynamics [13]. It incorporates individual elements in a simplest 
possible fashion as two-state stochastic oscillators with one quiescent state and one 
excited state, and random transitions between these two states which are influenced by 
the mutual coupling among the network elements. The model has been introduced in the 
deterministic limit of huge many coupled elements in complete neglect of the intrinsic 
mesoscopic noise and became immensely popular with the years [13], being used e.g. to 
describe neuronal oscillations in visual cortex within a mean-field approximation pun. 
Recently, the previously neglected mesoscopic noise effects were incorporated in this 
model for a finite-size network [151116]. Such a noise has been shown to be very 
important, in particular, for the occurrence of the critical avalanche dynamics na 
absent in the deterministic Wilson-Cowan model and also for the emergence of oscillatory 
noisy dynamics [16]. At the first glance, such a noisy dynamics can look like a chaotic 
deterministic one. Deterministic chaos influenced by the noise can also be a natural 
feature of a higher-dimensional dynamics, beyond the original two-dimensional Wilson- 
Cowan mean field model. Indeed, deterministic chaos has been found in the brain 
dynamics some time ago mm- However, it cannot be described within the memoryless 
Wilson-Cowan model because the minimal dimension for chaos is three ra. 

Stochastic mesoscopic noise effects due to a finite number of elements in finite size 
systems attract substantial attention over several decades, especially with respect to 
chemical reactions on the mesoscale PEGS], being especially pertinent to the physico¬ 
chemical processes in living cells [21]. In particular, such intrinsic noise can cause and 
optimize spontaneous spiking (coherence resonance [221123]) in the excitable clusters 
of ionic channels in cell membranes, which are globally coupled through the common 
membrane potential, or the response of such systems to periodic external signals 
(stochastic resonance [24]) within a stochastic Hodgkin-Huxley model [25]. Finite- 
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a) b) 



Figure 1. (a) Two-state model of neuron with one quiescent and one excited state, 
(b) Three-state model with one inactivated refractory state, (c) Extension of (b) 
incorporating multi-stage delay. 


size networks of globally coupled bistable stochastic oscillators were also considered 
without relation to the Wilson-Cowan model [261127] . including non-Markovian memory 

effects [231ESH3I!. 

In this paper, we consider a class of higher-dimensional generalizations of the 
stochastic Wilson-Cowan model aimed to incorporate non-Markovian memory effects 
in the dynamics of individual neurons. Such effects are caused by the existence of 
a refractory period or inactivated state from which the neuron cannot be excited 
immediately. First, we discuss a general class of such models. Then, we apply 
the simplest two-state non-Markovian model of this class, embedded as a three-state 
Markovian model with one inactivated state, to study a critical avalanche dynamics 
in a balanced network of the excitatory and inhibitory neurons within a mean- 
field approximation. Here, we restrict ourselves to the simplest example of a fully 
connected network with all-to-all coupling of its elements. In particular, we derive 
the power law exponents characterizing the critical self-organized dynamics of the 
network from the precise numerical simulations done with the dynamic Monte Carlo, 
or Gillespie algorithm. We also compare these results with similar results obtained 
within approximate stochastic Langevin dynamics, or, equivalently, within a diffusional 
approximation to the discrete state dynamics. Here, we reveal a profound difference. 
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2.1. Stochastic models of single neurons 


Let us depart from the Markovian model of a neuron possessing one activated or excited 
state “a”, and a quiescent state “q”, see Fig. [U a. The excitation of the neuron occurs 
with the rate f3f, where /3 is a rate constant and / is a dimensionless transfer function 
which depends on the states of the other neurons, and will be discussed below. Let 
us assume for a while that / is not explicitly time dependent. From the point of 
view of the theory of continuous time random walks (CTRWs) or renewal processes 
[32j[33], such a two state neuron can be completely characterized by the residence time 
distributions (RTDs) in its two states, ?/> a (f), an d ip q (t), correspondingly (assuming 
that no correlations between the residence time intervals is present - the renewal or 
semi-Markovian assumption). RTDs define completely the trajectory realizations of 
such a renewal process. In the Markovian case, the RTDs are strictly exponential, 
ifait) = aexp(— at), and f’qit) = z/exp(— vt), where we denoted v = /3f. Then, such a 
trajectory description corresponds to the Markovian balance or master equations for the 
probabilities to populate the states “a” and “q”, u(t) and q(t), correspondingly. Due to 
the probability conservation, u(t) + q(t) = 1, 

u{t) = — au(t ) + /3f[ 1 — u(t)] . (1) 


The memory effect due to a delay of a new excitation event after the neuron comes 
into the quiescent state, or the existence of some refractory period r^, can be captured 
within the trajectory description by a non-exponential RTD f> q (t). This transforms the 
corresponding master equation into a generalized master equation (GME) with memory, 
where the term /3fq(t ) is replaced by f* K(t — t')q(t’)dt ', with a memory kernel K(t). 
Here, t 0 is the starting time, f 0 = 0, if not a different one is explicitely stated. Hence, 
Eq. (jT|) is replaced by 


u(t ) = —au(f) + 


K{t — t')[l — u(t')\dt' . 


( 2 ) 


Jo 

In the CTRW theory it is well-known how the memory kernel K(t) and the residence 
time distribution 'ijjq(t) are related [33li3Tj (see also Appendix of [35]). Namely, their 
Laplace-transforms [denoted as F(s ) = f Q exp(—st)F(t)dt, for any function F(t)] are 
related as 


K(s) = 


S'fq(s) 


(3) 


1 - lfjg(s) 

In neurosciences, a delayed exponential, or delayed Poissonian model is popular 
It is featured by the absolute refractory period r^, i.e. il> q (t) = 0, for 0 < t < r^, 
and the exponential distribution, f>q(t) = z/exp[— v(t — r^)], for t > Td , see in Fig. 
[2] This model corresponds to 'ipq(s) = exp (—Tds)u/(s + is), and the memory kernel 
K(s) = vsj\{y + s) exp(rds) — v\. The numerical inverse Laplace transform of this 
memory kernel is depicted in the inset of Fig. [2l b. Notice that it does not correspond 





Figure 2. (a) Residence time distribution in the quiescent state, (b) and the 

corresponding memory kernel for the delayed exponential distribution and distributions 
corresponding to one, n = 1, and two, n = 2, inactivated states in Fig. [U b, and Fig. 
HI c, with n = 2 , correspondingly. Inset in part (b) shows also the cases n = 100, 
n = 1000, and n —> oo (delayed exponential). Time t is in the units of r d = I/ 7 , 
and v = 0.5. Numerical results in the inset were obtained by numerical inversion of 
the corresponding Laplace-transform using the Gaver-Stehfest method with arbitrary 
precision (35}[|D] to arrive at convergent results. 


to the memory kernel K(t) = u r d(t — T d ), which would correspond to the master equation 
with the deterministic delay [37] 

u(t) = -au(t) +Pfr[l - u{t - T d )] . (4) 


However, this memory kernel is strongly peaked at t = r d , and can thus be approximated, 
with v r = liin s _j.o K(s) — v /{1 + vr d ), which is the inverse mean time of the delayed 
Poissonian distribution ij> q (t ). I 11 the corresponding Markovian approximation, Eq. (JT]) 
is restored with a renormalized transfer function, 

/ 


fr = 


(5) 


1 + /3r d f 

This is the simplest way to account for the delay effects. Obviously, any delay should 
suppress excitability within this approximation, because f r < f. However, suppression 
of the excitability of the inhibitory neurons may enhance the excitability of the whole 
network consisting of both excitatory and inhibitory neurons. Hence, possible effects are 
generally nontrivial even in this approximation. Moreover, Eq. (J5]) makes it immediately 
clear that the delay effects are generally expected to be very substantial for /3T d > 1 . 


2.1.1. The simplest non-Markovian model and its Markovian embedding. It is well- 
known that in many cases non-Markovian CTRW dynamics can be embedded as some 
Markovian dynamics in a higher-dimensional, possibly infinite dimensional space [38]. 
Given a non-trivial form of the memory kernel for the delayed exponential distribution 
of the quiescent times, we can ask the question: What is the simplest non-Markovian 
model and the corresponding Markovian embedding to account for the memory effects? 
From the point of view of GME, it is K(t) = Kexp(— rt), i.e. an exponentially 
decaying memory kernel. The corresponding memory kernel with k = 1 / 7 , and 
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r = v + 7 corresponds to i/j, 7 {t) , which is the time convolution of two exponential 
distributions, i/jq°\t) = uexp(-ut), and ^(t) = 7 exp(— 7 1). It corresponds to a 
compound state “q = q U i” in Fig. [U b. Indeed, the Laplace transform of the 
corresponding compound distribution is just the product of the Laplace-transforms 
of two exponential distributions, i.e., 'ipq(s) = 1 / 7 /[(s + u)(s + 7 )]. By Eq. (J3]) it 
corresponds precisely to the stated exponential memory kernel. The corresponding 
ip q (t) = u r y[exp(—ut) — exp(— 7 t)\/{v — 7 ) has a maximum at the most probable time 
interval f max = ln(u /'y) / (u — 7 ), see Fig. [2l a, reflecting the most probable stochastic 
time delay. This simplest non-Markovian model with memory allows, however, for a 
very simple Markovian embedding by introduction of an intermediated refractory state 
“i” shown in Fig. [D b, with the population probability x(t) and the exponential RTD 
given above. It has the mean refractory time 77 := (r) = ./ 0 °° V ; *( r )^ r = 1/7, and the 
relative standard deviation, or the coefficient of variation, Cy := \J (r 2 ) — (r) 2 / (r) = 1 . 
Using the conservation law, u + q + x = 1, the corresponding master equations can be 
written either as 

u — — au + v{l — u — x), 

x = au — 'yx, (6) 

or as 

u — — au + uq, (7a) 

q = - (7 + u)q + 7(1 - u) . (7b) 

From (17b|) follows 

q(t) = e-^ +v)t q( 0 ) + 7 f - u{t')]dt' . ( 8 ) 

Jo 

After substitution of this equation into (I7aj) one obtains indeed Eq. (J 2 J) with the 
discussed exponential memory kernel provided that q{ 0) = 0. The latter condition 
is natural because every sojourn in the compound quiescent state “q” starts from the 
substate “i” (resetting memory of this neuron to zero), and q(t) +x(t) is the probability 
of the compound quiescent state within the two-state non-Markovian reduction of the 
three-state Markovian problem. Here, one can also see the origin of a profound problem 
with the description of the whole network dynamics of interacting non-Markovian 
renewal elements as a hyper-dimensional renewal process. Obviously, the behavior of the 
whole network cannot be considered as a renewal process, because after each and every 
de-excitation event only one element is reset. Then it starts with zero memory, while 
all others keep their memory until they are reset. Hence, any Gillespie type simulation 
of the whole network of interacting non-Markovian elements must account for the “age” 
of each network element separately. Markovian embedding allows to circumvent this 
problem and dramatically accelerate simulations within the mean-held approximation, 
see below. 

The considered three-state Markovian cyclic model presents one of the fundamental 
kinetic models in biophysics. It provides, in particular, a paradigm for non-equilibrium 
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steady state cycling. For example, the cyclic kinetics of an enzyme E , which binds 
a substrate molecule S, converts it to a product P, and releases it afterwards can 
be represented as a three-state cycle, E —> ES —$■ EP —* E. This model was used 
e.g. in Ref. [29] for an excitable unit. Furthermore, three-state non-Markovian models 
can be used with a non-exponential distribution ipi{t). For example, if to use the 
deterministically delayed ^(t) = S(t — r^), and exponential ijj q [t) within the three- 
state cyclic model, then one obtains the delayed exponential distribution within the 
two-state reduced model, which was discussed above. In addition, i/} a {t) can also be 
non-exponential. For ^ a (f) = 5{t — 1/a) in the excited state of the three-state non- 
Markovian model, one obtains the model used in Refs. [281129] . 


2.1.2. Markovian embedding with many substates. One can also introduce many 
delayed substates as shown in Fig. [I] c. Within the three-state non-Markovian model 
this can be considered as having one delayed state “i” characterized by the special 
Erlangian distribution [32], ipi(t) = n^{n^t) n ~ l exp(— njt)/{n — 1)!, with the Laplace- 
transform = 1/(1 + s/{n^()) n reflecting the corresponding multiple convolution. 

Such a non-Markovian three-state model has been considered in Bm. and a non- 
Markovian two-state model with the Erlangian distribution of the quiescent times has 
been studied in [31]. The compound quiescent state corresponding to the model in [30] is 
characterized by i^ q {s) = u /[(1 + s / {n^)) n {y + s)]. The mean delay time is the same 77 = 
l/y for any n, and the coefficient of variation becomes ever smaller with increasing n, 
Cy = 1 /y/n, i.e. the distribution of the refractory times becomes ever more sharpened. 
The Laplace-transformed memory kernel is K(s) = vs/[{ 1 + s/{n'y)) n {y + s) — v\. Some 
corresponding ip q (t) and K{t) are shown in Fig. [2J Already for n — 2, the memory 
kernel starts to show a peaked structure. Notice that in the limit n —?• 00 , the above 
delayed exponential (or Poissonian) model immediately follows with 77 = I/ 7 . For any 
n, the inverse mean time in the quiescent state is given by /3f r with f r in (]5]h Increasing 
n yields an ever better approximation for the delayed Poissonian model. However, it 
can be considered as a useful model in itself. The corresponding Markovian embedding 
master equation reads (with q excluded by the probability conservation law): 


u = — era + (3f 


Xi = au — 77 . 7 x 1 , 



Xj = n'.ix, - Xj), j = 2, ,.,n, 


(9) 


with, Xj(0) = 0, for j = 2, ..., 77 , initially. With a different initial condition, the 
corresponding GME obtained upon projection of the multi-dimensional dynamics onto 
the subspace of u and q variables will contain a dependence on this initial condition in 
the subspace of hidden Markovian variables. On the level of non-Markovian dynamics 
this can be accounted for by a different choice of the residence time distribution 'i/j q 0 ' > (t) 
for the first sojourn in the quiescent state. It depends on how long this state has been 
populated before the dynamics started [32]. The GME (J2J) , (j3J) corresponds to the 
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particular choice, ipq°\t ) = fi’qif)- 

We mention in passing also that it is straightforward to consider a power-law 
distributed delay, both within a semi-Markovian model and within an approximate 
finite-dimensional Markovian embedding. Also a stochastic model for bursting neurons 
can be introduced immediately. We shall not, however, consider these possibilities in 
the present work. 

2.2. Network of neurons within the mean field dynamics 

Following Wilson and Cowan, we consider a network of N e excitatory and N, t inhibitory 
neurons, with the probabilities of neurons to be in their excited states ufit) and vfit), 
correspondingly. The neuron k can influence the neuron l and possibly itself {k — l), 
by excitation, or inhibition with the coupling constants, w l fi e > 0, w l fi > 0 for the 
excitatory neuron k, and —w l fi < 0, —w\f < 0, for the inhibitory neuron k. The absolute 
value of the coupling constant reflects the synaptic strength. Each excitatory neuron l 
thus obtains an averaged input si = (1 /N l e ) Yk w ee u k — (1 /N\) fifip w l fiv p + h l e , and the 
inhibitory neuron p receives the input s p = (1/Nf) Yi wf e ui — (1 /Nf) Y2k w ii v k + h%, 
where N l e and N\, etc. is the number of inputs which the Z—tli neuron obtains from the 
excitatory and inhibitory neurons, correspondingly. The constants h l e and h\ serve 
to fix the spontaneous spiking rates, fiif(h l e ), fipffihfi), in the absence of coupling, 
w l fi e —> 0,w l £ —> 0,w l fi —$■ 0,w l £ —> 0. Coupling can either enhance, or suppress these 
rates. Phenomenologically, this is accounted for by the transfer function /(s), which we 
assume to be the same for all neurons. Some common biophysically motivated popular 
choices of the transfer function /(s) are 

/(s) = tanh(s)#(s), (10) 

where 9(s ) is the Heaviside step function, and 

f(s) = l/[l + exp(-s)] . (11) 

Both are bounded as 0 < / < 1. Evidently, this is a very rich model even for the simplest 
two-state model of neurons, in the absence of memory effects. The simplest further 
approximation to describe the collective dynamics of neurons is to invoke the mean 
held approximation [14]. It is equivalent to assuming that all the coupling constants 
like w l fi e , etc., thresholds h l e , etc., and rates fii, on , are equal within a subpopulation, 
w l fi e = w ee , h l e = h e , /3 l e = fi e , or fi\ = fii, etc. Furthermore, one can introduce the 
occupation numbers of the excited neurons in each population, u{t) = (1 /N e ) Y^i=i 
and v(t) = (1 /Nf) anc i consider the dynamics of these variables. They present 

the fractions of neurons which are excited. 

We restrict our treatment in the rest of this paper to the simplest two state 
non-Markovian model within the three state Markovian embedding and introduce the 
occupation numbers of neurons, x(t ) = (1 /N e ) an d y{t) — (1 /Nf) 

in the corresponding delayed states. Then, in the deterministic limit N e ,Ni —> oo, one 
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obtains a 4-dimensional nonlinear dynamics, 

u = — a e u + (3 e f (w ee u - w ei v + h e )(l - u - x), 
x = a e u — 7 e x, 

v = - aiV + Pif(w ie u - wav + hi)( 1-v - y ), 
y = Oi iV - 7 iV • ( 12 ) 

Notice that unlike the original two-dimensional mean-held Wilson-Cowan dynamics 
in the deterministic limit, the considered 4-dimensional dynamics can in principle be 
chaotic, for some parameters (which remains an open question). Dynamical chaos might 
emerge already when only one sort of neurons, e.g. inhibitory neurons, exhibits delayed 
dynamics, since the minimal dimension for nonlinear chaotic dynamics is three. Then, 
in the macroscopic deterministic limit, 

u = - a e u + f3 e f (w ee u - w ei v + h e ){ 1 - u), 
v = — otiV + Pif(w ie u - Wav + hi)( 1-v- y), 

y = a iV - 7 iV . (13) 

However, we shall not investigate the possibility of a deterministic chaos emerging due 
to a delay within the minimal extensions of the Wilson-Cowan model in the present 
work, but rather focus on the mesoscopic intrinsic noise effects caused by finite N e and 
i\. Then, the occupational numbers are random variables (at any fixed time t). 


2.2.1. Langevin dynamics. For a very large number of neurons, one can account for the 
mesoscopic noise effect within the Langevin dynamics, or the diffusional approximation 
of the discrete state birth-and-death process describing the evolution of the network. 
This procedure is standard, by analogy with the stochastic theory of chemical reactions 
m- Since we have only direct “reactions” like q —> a, a —> i, i —» q, for two type of 
neurons, one must introduce six variables and six independent zero-mean white Gaussian 
noise sources £,(f), (£i(t)£j(t / )) = 5ij5(t — t'). Stochastic dynamics is, however, effectively 
4-dimensional because of two probability conservation laws, which allow to exclude two 
variables out of six: 


u = — a e u + /3 e f(w ee u — w e iV + h e )(l — u — x) 

--^a^U^iit) + — 7 = \J fief (u! ee U - W ei V + h e ){ 1 ~U~ x)£ 2 (t), 
v-^e 

1 ___ 1 




x = a e u — 7 e x + —=y/a^j^(t) -7=Vt^ f 3 (*), 

YiV e y iV e 

v = - aiV + [dif (w ie u - WaV + hi)( 1 -v- y) 

faivU(t) + — r =\/Pif(w ie u - wav + hi)(l - v - y)£ 5 (t), 

\ i 


y/Nl 


y = a iV - liV + -^v^T^ 4 (t) - . 


(14) 


In the limit N e ,Ni —y oo, the deterministic description in Eq. (1T2]) is restored. The noise 
is multiplicative and the Langevin equations must be Ito-interpreted, as it is always the 
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Table 1. Transitions and rates 


i 

Transition 

Rate ri 

1 

(N, M, Ni, Mi) -> (N — l,M,Ni + l,Mi) 

r i = Na e 

2 

(N, M, Ni, Mi) ->■ (N,M-1, Ni,Mi + l) 

r 2 = Mai 

3 

{N, M, Ni,Mi) -+ (N + 1, M, Ni, Mi) 

r 3 = (N e -N- Ni)p e f( Wee N/N e - 
w ei M/Ni + h e ) 

4 

(N, M, N u Mi) —► ( N,M+ 1, Ni, Mi) 

U = (Ni-M- Mi)/3if(w ie N/N e - 

wuM/Ni + hi) 

5 

(N, M, Ni,Mi) -> (N, M, W - 1, Mi) 

r 5 = 7 e Ni 

6 

(N, M, N h Mi) (N, M, N h Mi- 1) 

r 6 = 7 iMi 


case if the Langevin dynamics results from the standard diffusional approximation to a 
birth-and-death process, or chemical master equation [20]. Notice that such a Langevin 
stochastic description can become problematic, if any of the variables u,v,x,y becomes 
temporally zero or one. Even if some of the noise terms do vanish at the boundaries, 
where there corresponding rates vanish, the others do not, when a particular boundary 
is hit. Hence, the occupational numbers can in principle become temporally negative, 
or larger than one. This unphysical feature is produced by the standard diffusional 
approximation. However, this problem can be fixed in the numerical simulations by 
introduction of the corresponding reflecting boundaries and taking sufficiently small 
integration time steps, as done e.g. in Ref. [251I4T] for stochastic Hodgkin-Huxley 
equations. 

2.3. Exact stochastic simulations. 

Within the mean-field approximation of Markovian dynamics, it suffices to count the 
numbers of neurons in the corresponding activated, N and M, and refractory, N\ 
and Mi states. Then, we are dealing with a random walk on a 4-dimensional lattice 
(N, M, N\, Mi) with the discrete variables N, and N\ taking values in the range from 
zero to N e , and the variables M and Mi in the range from zero to W, so that also 
0 < N + Ni < N e and 0 < M + Mi < Ni. From the site (N, M, Ni,Mi) six 
different transitions are possible. They are enlisted in Table CD with the corresponding 
transition rates. The master equation governing this birth-and-death process can be 
readily written. However, it is bulky and not very insightful. For this reason, it is not 
presented here. The corresponding stochastic process can be easily simulated with the 
dynamical Monte Carlo or Gillespie algorithm ra, which is exact. Namely, on each 
step one draws two random numbers. The first one, r, is drawn from the exponential 
distribution characterized by the total rate = Xu =1 W It gives a random time interval 
at which the network state is updated. Given a uniformly distributed random variable 
Cl; 0 < Ci < 1, r = (l/r s ) ln(l/Ci). Then, one of the transitions in Table [l] is chosen 
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in accordance with its probability p r = ?y/r£. For this, one generates a uniformly 
distributed random variable ( 2 bounded as 0 < £2 < 1- If 0 < £2 < Pi> then the first 
transition is chosen. If pi < ( 2 < Pi + P 2 , then the second transition is chosen, etc., i.e. 
in accordance with the length of the corresponding interval p t , Y^!i= 1 Pi = 1- 

Notice that an attempt to generalize this scheme towards a non-Markovian renewal 
walk on a 4-dimensional lattice to account for the memory in the inactivated state 
is logically inconsistent because in such a case accomplishing each step would mean 
reset, or renewal of all neurons, and not the only one which actually makes transition. 
However, each non-Markovian element has its individual memory. In a direct simulation 
of the network of non-Markovian elements one must therefore consider them individually, 
even within the mean-field approximation. Then, one has to consider CTRW on a 
hyper-dimensional lattice of huge dimensionality, which will dramatically slow down 
simulations imposing computational restrictions on the maximal size of the network. 
Of course, beyond the mean field approximation one must also simulate each element 
separately. Here, a direct semi-Markovian approach can be preferred. In this work, 
we restrict ourselves to the mean-held dynamics within the Markovian embedding 
framework, which allows for exact simulations of very large networks within a reasonable 
computational time. 

3. Results and Discussion 

3.1. Oscillatory dynamics of neuronal network 

First, we test our stochastic simulations done with XPPAUT ia against nonlinear 
deterministic dynamics for a very large network size with N e = 10 6 and = 10 6 . 
For this, departing from the parameter set in Ref. JT 6 ] (the case of without delay) we 
use a set of parameters, where an oscillatory dynamics emerges: a e = 0 . 1 , cc* = 0 . 2 , 
fie I 7 fii 2, y e 10, y i 10, Wee 32, W e i 32, Wie 28, Wa 2, he 3.8, 
hi = —9, and the transfer function in Eq. (TTll) . Time is in milliseconds and the rate 
constants are in inverse milliseconds. The difference is barely detectable in Fig. [3j a, b, 
where we present the results of stochastic simulations done both with the exact Gillespie 
algorithm and within the approximate Langevin dynamics. However, stochastic effects 
become immediately seen in Fig. [3j c, d, where we reduced the number of neurons to 
N e = 10 3 and iV* = 10 3 . We also compare in Fig. [3j a, b, the results for the considered 
dynamics and its two-variable Markovian approximation given by the standard Wilson- 
Cowan model in which, however, the transfer functions are renormalized in accordance 
with Eq. ([5]), where the parameter fird is replaced by fi e / y e = 0.1 and fii/^i = 0.2, 
correspondingly. The deviations are visible, but small. However, the differences become 
very pronounced for small y e = y* = 0.1 corresponding to the mean refractory period 
Td = 1/y e ,i of 10 msec. Then, the Markovian approximation fails completely, see in Fig. 
|4] especially in part (b), revealing that neither the form of the oscillations, not their 
period are reproduced even approximately. Especially remarkable is that contrary to 
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Figure 3. (a) Limiting cycle in the u — v plane for deterministic dynamics and for 
stochastic dynamics, u = N/N e , v = M/Ni, with N e = Ni = 10 6 , (b) Time-dependence 
of the u variable for (a), (c) Limiting cycle in the u—v plane for deterministic dynamics 
and for stochastic dynamics with N e = Ni = 10 3 , and (d) the corresponding time- 
dependence of the u variable. Langevin simulations are done with the stochastic Euler 
algorithm using time step St = 0.05 in (a) and (b), and St = 10 -5 in (c) and (d). 

intuition the increase of the refractory period of a single neuron does not increase the 
period of oscillations, as the Markovian approximation predicts, but rather makes it 
smaller - the tendency is opposite! Therefore, non-Markovian memory effects generally 
do matter and one should take such effects seriously into account. With a small further 
decrease of 7 e , 7 i to 7 e = 7 i ~ 0.08873 with 77 & 11.270 the oscillations are terminated 
by a supercritical Hopf bifurcation (not shown). Interestingly, Markovian approximation 
also predicts such a termination, but at a slightly larger critical value 7 e = 7 i ~ 0.0987 
with critical 77 ~ 10.132. This makes clear that the phase transitions between the 
quiescent network and the network undergoing synchronized oscillations are possible 
with respect to the length of the refractory period used as a control parameter. 

3.2. Noise-induced critical avalanche dynamics 

In the remainder, we investigate the influence of memory effects on the avalanche 
dynamics. As it has been shown in [15], in order to have avalanche dynamics the 
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Figure 4. Deterministic two-state dynamics with memory and its Markovian 
approximation for 7 e = 7* = 0.1. Other parameters are the same as in Fig. [3] 
The deviation in (b) indicates that the Markovian approximation fails completely. 


excitatory and inhibitory processes should be nearly balanced, and the network should 
have a so-called feedforward structure. Then, one can achieve a sort of self-organized 
critical (SOC) state |T| [ 8 j sustained due to intrinsic mesoscopic fluctuations. Very 
different from other SOC models, here fluctuations play a major role and in the 
deterministic limit avalanches disappear, i.e. they are of mesoscopic nature. The 
nullclines of 2 d deterministic dynamics in the absence of memory effects, u = 0 and 
v — 0, should cross at a very small angle in the u — v plane, so that fluctuations can 
produce large amplitude outbursts of the u and v variables moving synchronously but 
randomly, i.e. the subpopulations of excitatory and inhibitory neurons are synchronized 
exhibiting stochastic dynamics at the same time [15]. In the same spirit, we choose 

a e = Oii = 1, & = /?» = 5, w ee = W ie = w e = 30, W e i = Wu = w t = 29.9, so 

that w e — Wi <C w e + Wi, and overall excitation slightly dominates over inhibition. 
Furthermore, we choose h e = hi = 0.001, and the transfer function in Eq. (fTOjh as in 
Ref. [15]. The rates 7 e , 7 \ and the number of neurons were varied. Large 7 e = 7 * = 10 
corresponds to a small refractory time of 0.1 msec, whereas 7 e = 7 * = 0.1 corresponds to 
a profound delay with 77 = 10 msec, so that the individual spiking rate of neurons cannot 
exceed 100 Hz being limited by the refractory period. Typical avalanche dynamics is 
shown in Fig. [5] for L(t) = N(t ) + M(t ) with N e = Ni = 10 3 , and (a) 7 e = 7 \ — 10, 

(b) 7 e = 7 i = 1, (c) 7 e = 7 i = 0.1. Furthermore, in Fig. [5] d, e, f, we show the 

influence of an increasing number of neurons on the avalanche dynamics. The following 
tendencies are clear. First, the increase of refractory period reduces the maximal amount 
of neurons involved in spiking, from about 76% in Fig. [5] a to 58.5% in Fig. [5j c. Such 
a tendency is already expected from the renormalization of the transfer function in the 
Markovian approximation, cf . Eq. (J5]) . However, this tendency is in fact much weaker 
since / 3 t 7 = 50 in the part (c), and the renormalization argumentation would predict 
almost complete suppression of avalanches for such a delay. Even more astonishing is 
that avalanches still did not vanish even for a very large = N e = 10 6 , see Fig. [5j f. 
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Figure 5. Avalanches for N e = Ni = 10 3 and (a) 7 e = 7 i = 10, (b) y e = 7 » = 1, 
(c) 7 e = 7 i = 0.1. In (d), (e), and (f), 7 e = 7 i = 10, however, the network size 
is increased to (d) N e = Ni = 10 4 , (e) N e = Ni = 10 5 , and (f) N e = Ni = 10 6 . 
For the fixed N e = Ni = 10 3 , the maximal L(t) in (a) is 1513, i.e. about 76% of 
maximally possible. With the increase of refractory time it diminishes to 1322 (about 
66 %) in (b), and to 1169 (about 58.5%) in (c). For a fixed refractory time, but with the 
increase of the network size the maximal number of active neurons is 12791 (~64%) 
in (d), 61923 (~31%) in (e), and 223234 (~11%) in (f). With the increase of network 
size, the relative size of avalanches decreases. Notice also that the minimal number 
of active neurons is L m ; n = 3 in (e), and L m j n = 112 in (f). This must be taken 
into account when one defines avalanches in large networks. Otherwise, one can come 
to incorrect conclusion that the avalanches cease to exist, which is manifestly refuted 
in (f) for a very large number of neurons which seems to be macroscopically large, 
and nevertheless fluctuations are still very important, even though they do vanish 
in the strict limit N e ,Ni oo. Experimentally, one also defines the start and end 
of an avalanche by crossing a threshold of basal activity upwards, and downwards, 
correspondingly. Simulations are done with the Gillespie algorithm. 


This is very different from the oscillatory dynamics of a network of the same size, cf. 
Fig. |3l a, and b, which is practically deterministic. Of course, with increasing network 
size, the relative amplitude of avalanches becomes ever smaller, and there also emerges 
a minimal number of neurons excited, i.e. the network activity never goes down to zero. 
However, this is also so in the real neuronal dynamics. Such a dominance of mesoscopic 
fluctuations in a system of millions elements with a special (feedforward) structure of 
coupling is really astonishing. This is a feature of some critical state, as we know from 
statistical physics. 

To statistically characterize the avalanche size distribution and their duration we 
proceed in accordance with the procedure outlined in Ref. [15]. It reflects, in part, 
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Figure 6. Survival probabilities for the avalanche size (a), duration (b), and peak 
(c) distributions obtained from the Gillespie simulations. The case without delay. 
Circle symbols correspond to numerical results and lines to the corresponding fits 
with parameters displayed on the plot. Blue line corresponds to the power law fit. 
N e = Ni = 10 3 . Other parameter are given in the text. 


also the experimental procedure cm- Namely, we first discretize the time series with a 
time bin of the size (At), which corresponds to the averaged interspike time distance 
in a particular simulation. Then, an avalanche is defined by its start, when the spiking 
activity crosses some threshold level L t h r , and its end, when the network activity drops 
to (L t h r = 0) or below L t ^ r > 0 after some time, which defines the avalanche duration. 
The size is defined as the sum of the number of neurons active in each time bin during 
the avalanche. It is also defined experimentally in such a way. In essence, the size S of 
an avalanche is the integral of the network activity in Fig. [5] over the time during each 
avalanche period divided by the time bin width. Of course, as also in experiments the 
critical exponents discussed below depend both on the time bin width and on the basal 
level of neuronal activity L t h r . However, this dependence is weak for a truly critical 
dynamics. By doing statistical analysis, we first find the survival probability F(S), or, 
equivalently, the cumulative probability 1 — F(S) from the numerical data. Then, the 
distribution density follows as p(S) = —dF(S)/dS. 
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Let us start from the case without any time delay, 77 = 0. The survival probability 
for the avalanche size distribution F(S) is shown in Fig. [Gja, for the time bin 
(At) = 0.00643597 and L t h r = 0. It shows three characteristic features: (1) an initial 
Weibull distribution, F(S ) = exp(— [(S — l) 01 /^ 1 ]); with cp « 0.71 and 5j « 670; (2) 
an intermediate power law F(S) oc S a2 , with 02 ~ —0.165, and (3) an exponential tail 
F(S) = pi exp(—S/Si) with pi ~ 0.204 and S 2 = 6.97 x 10 5 . The size distribution p(S) 
is, therefore, initially approximately a power law with negative exponent a\ — 1 ~ —0.29, 
followed by a power law with negative exponent a = a 2 — 1 ~ —1.165. The latter 
one extends over approximately two size decades and ends with an exponential tail 
characterized by a cutoff size, Si. The corresponding survival probability F(T) for the 
avalanche durations T is shown in Fig. [ 6 ] b. It can be well fitted by a sum of two 
exponentials, 

F(T) = pq exp(—T/T 0 ) + (1 -p 0 )exp(-T/Ti) . (15) 

However, it also seems to display an intermediate power law over about one time decade, 
from 1 to 10 ms, with the power exponent 61 ~ —0.523, and the cutoff time Tj ~ 10.93 
ms. Hence, the probability distribution p(T) = —dF(T)/dT also appears to reflect an 
intermediate power law p(T) oc T b with b — b\ — 1 ~ —1.523. Interestingly, the duration 
of avalanches in experiments with organotypic cortical neuronal systems has a similar 
cutoff time of about 10-20 ms, with a maximal avalanche duration of about 40-80 ms, 
which is restricted by the period of 7 —oscillations ms- The intermediate power law 
also extends over about one time decade in the experiments. However, the experimental 
power law exponent is different, b exp fx —2. It should be mentioned in this respect that 
the time bin in the experiments is also very different, At ~ 1 — 4 ms. One electrode 
measures in experiments a contribution of many neurons. In fact, coarse graining over 
some unknown A L should be done. The experimental size exponent a is also different, 
otexp ~ —1.5. It is not, however, the goal of this paper to provide a model fully consistent 
with the experimental observations, which are subject of ongoing research work and some 
controversy in the literature [43]. In this respect, a bi-exponential dependence can be 
perceived as a power law over one intermediate time decade, as our fit also shows. 

As an additional characteristics of the avalanches size, one can also consider the 
maximal number of neurons activated at once during an avalanche, or the avalanche 
peak with a distribution density p max (S) = —dF max (S)/dS. The corresponding survival 
probability, F max (S ), also exhibits a power law, F max (S ) oc S a ' 2 , with a' 2 ~ —0.175, in 
Fig. [ 6 ] c. Hence, p max (S) oc S a \ with a' = a 2 — 1 ~ —1.175, which only slightly 
differs from a ~ —1.165 meaning that the avalanche size is roughly proportional to 
its peak. However, the cutoff of F max (S ) is super-exponentially sharp, because the 
maximal number of neurons involved in an avalanche at the same time is restricted 
by the total number of neurons in the network. Furthermore, F max reveals a very large 
portion of avalanches whose peak does not exceed 10 , which explains the initial stretched 
exponential dependence in Fig. [ 6 ] a. Strictly speaking, this part of the size distribution 
(with L t h r = 0 ) reflects a background or basal noise, where neurons practically do not 
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Figure 7. Survival probabilities for the avalanche size (a), duration (b), and peak (c) 
distributions for the network with a time delay, 7i = 7 e = 10. Other parameters are 
the same as in Fig. [G] The cutoff size of avalanches Si becomes slightly smaller than 
in Fig. [Gj and the characteristic power law exponents are also slightly changed. 


interact with each other, and there are no avalanches of spontaneously increased activity, 
which are characterized by a power law distribution. 

Next, we like to clarify how robust these features are for networks with a time 
delay. For this, we study the influence of the mean delay time by decreasing the rates 
7 e = 7 i from 10 through 1.0 to 0.1 in Figs. [7] 0 El respectively. The mean delay time 
increases, accordingly, from 0.1 through 1.0 to 10 ms. Even though the parameters 
of the distributions do change, these changes are not dramatical. In particular, the 
corresponding critical size exponent a changes from —1.165 (no delay), to —1.152, 
— 1.171 and —1.161, respectively. Accordingly, the critical time exponent b changes 
from —1.523 (no delay) to —1.505, —1.558, and —1.511, respectively. Such changes 
are not statistically significant, and one cannot detect any systematic tendency upon 
a variation of 77 . The point is that these exponents are also changed a bit, if we use 
e.g. 2(A t), or 3(A t) for the time bin (not shown). They also depend on the threshold 
L t hr- In this respect, if to change L t h r from 0 to 10, the initial stretched exponential 
part of the size distribution practically disappears. However, there appears an initial 
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Figure 8. The time delay is increased further with respect to the one in Fig. 0 Here, 
7 * = 7 e = 1. The other parameters remain the same. The cutoff size of avalanches, 
Si = 6.33 x 10 5 , is now visibly smaller than one without delay, Si = 6.97 x 10 5 in 
Fig. [6] The cutoff time Ti = 12.46 is increased with respect to Ti = 10.93 in Figs. 
10 i.e. the avalanches last longer. The power law exponents here deviate slightly in 
the opposite direction from the one in Fig. 0 They become closer to the case without 
delay in Fig. 1 This indicates that the time delay does not affect significantly the 
power law exponents. 


power law instead, see in Fig. [TUI Remarkably, the intermediate power law exponent 
remains rather robust. It is changed from a = —1.171 in Fig. [71a to a = —1.144 in Fig. 
COl a. This is a small variation. Notice, however, that the results in Fig. [TQ1 b in fact 
reject the hypothesis that there is an intermediate power law in the time distribution of 
avalanches. First, the power law region changes from larger to smaller times, and also 
(more important!) the corresponding time exponent changes from b\ = —0.558 in Fig. 
d b to b\ = —0.221 in Fig. [TUI b. Clearly, such a strong influence of the choice of L t hr 
on the “power law” exponent b makes it clear that this is not a power law. In fact, the 
time distribution is clearly biexponential. 

Though plausible until this point, it remains, however, strictly speaking, still not 
quite clear if a is indeed a critical exponent. If true, the extension of the power law 
domain of the whole size-distribution and the cutoff size S\ should increase with the 
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Figure 9. Influence of a further increase of the time delay by an order of magnitude, 
Td = 10, on the distributions depicted in Figs. [3 [8] Here, 7i = 7 e = 0.1, and the 
other parameters are not changed. Si drops further to Si = 6.07 x 10 5 , and T± slightly 
increases to T\ — 12.65. The power law exponents exhibit, however, merely some 
fluctuations without any systematic trend in Figs. EK] 


system size accordingly. Indeed, if we increase the system size tenfold keeping the other 
parameters the same as in Fig. [9l the power law domain in the size distributions also 
increases by an order of size magnitude, see in Fig. [TU a. The time cutoff Tj also 
increases in Fig. [TU b, i.e. the avalanches become longer. Also with decreasing the 
system size tenfold the power law domain shrinks accordingly in size, see Fig. [T2l 
a, and the avalanches become essentially shorter, as indicated by the decreased cutoff 
time Ti in Fig. [T21 b. Such scaling dependencies on the system size are typical in 
experiments. From this we can conclude that the size exponent a is indeed a critical 
exponent. However, within the considered model the avalanches do gradually vanish 
with an increase of the system size. Therefore, the adjective “critical” should be used 
also with respect to the exponent a with some reservations. We consider a rather 
atypical SOC model, even though the exponent a is by chance close to that of the 
sandpile model m- It should also be noticed that the initial distributions of the sizes 
and times and the tail functional dependencies can be sensitive to both the system 
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Figure 10. Influence of the choice of the detection threshold L t hr on the distributions 
of avalanche sizes (a) and time durations (b). Here, L t hr = 10 is used for the data 
analysis instead of Lthr = 0 in Fig. [3 for the same data. Noticeably, the initial 
stretched exponential part of the size distribution in Fig. [3 a disappears. Instead, 
there appears initially another power law dependence. The intermediate power law 
exponent 02 is, however, pretty robust, <X 2 = —0.144 here versus 02 = —0.171 in Fig. 
[3 In contrast with this, the intermediate power law exponent in the time distribution is 
changed dramatically from b\ = —0.558 in Fig. [3 to &i = —0.221. This fact disproves 
the hypothesis of an intermediate power law in the time distribution. It is clearly 
bi-exponential, Eq. ©• 



Figure 11. Influence of the increased network size on the distributions of the avalanche 
sizes (a) and the time durations (b). Here, Ni = N e = 10 4 vs. Ni = N e = 10 3 in Fig. 
El The other parameters are the same. The power law regime in the size distribution 
extends by an order of magnitude, with the cutoff size increased to Si = 6.54 x 10 6 , 
accordingly. The corresponding power law exponent varies insignificantly. The time 
cutoff T\ increases in (b) to T\ = 16.02 from T\ = 12.65 in Fig. [3 i.e. the avalanches 
last longer. 
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Figure 12. Influence of the decreased network size on the distributions of avalanche 
sizes (a) and time durations (b). Here, = N e = 10 2 vs. Ni = N e = 10 3 in Fig. [9] 
Other parameters are the same. The power law regime in the size distribution shrinks 
by an order of magnitude, with the cutoff form changed from the exponential in Fig. 
El a, to the Gaussian here. The intermediate power law exponent is not changed, 
however, strongly. The time cutoff Ti decreases in (b) to T\ = 4.79 from Tf = 12.65 
in Fig. [9] b, i.e. the avalanches become significantly shorter. 


size and the choice of the threshold L thr . For example, the size distribution exhibits a 
Gaussian tail in Fig. [T2l a, for a small system size. The intermediate power law in the 
size distribution is, however, rather robust, with a being in the range [—1.207, —1.144] 
for the data presented, with the average (a) ~ —1.164. 

3.2.1. Langevin dynamics of avalanches. Within the Langevin approximation of the 
discrete state dynamics, the avalanches look very similar to the ones depicted in Fig. 
[5j However, their statistics is very different. We performed the corresponding Langevin 
simulations for the same parameters as in Figs. [3 [TO] with the time step taken to 
be the same (At) = 0.00618608 as the time bin used to produce the results in Figs. 
m n We also used L t h r = 10 to analyze the data, as in Fig. |T0l The results 
shown in Fig. [13] reveal similar intermediate power laws both in the size and the 
time distributions yielding ol ~ —1.026, 6^ ~ —1.058. However, these results differ 
essentially from the results obtained within the exact dynamic Monte Carlo simulations, 
compare Fig. [13] with Fig. [TO] This indicates that the Gauss-Langevin or diffusional 
approximation of the genuine discrete state dynamics can deliver incorrect results for 
the fluctuation-induced avalanche dynamics. This fact makes any analytical theory 
for the numerical results presented in this work especially challenging. It is almost 
hopeless to develop such a theory for the discrete state avalanche dynamics within the 
studied model. Multi-dimensional birth-and-death processes are very difficult for any 
analytical treatment. Within the Langevin dynamics approximation, or the equivalent 
multi-dimensional Fokker-Planck equation description an analytical treatment is more 
feasible. However, such a theory will not help to understand the critical features of the 
discrete state dynamics, as our numerical results imply. 
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Figure 13. The results derived from the Langevin dynamics simulations for the size 
(a) and time (b) distributions, at the same parameters as in Fig. 1101 obtained with 
Lthi = 10. Notice the dramatical changes. 


4. Summary and Conclusions 

In this paper, we studied a generalization of the stochastic Wilson-Cowan model of 
neuronal network dynamics aimed to incorporate a refractory period delay on the 
level of individual elements. Considered as stochastic bistable elements such model 
neurons exhibit non-Markovian dynamics with memory, which can be characterized 
by a non-exponential residence time distribution in the resting state of the neuron 
(semi-Markovian description), or, alternatively, by the related memory kernel within a 
generalized master equation description. Such a non-Markovian description generally 
allows for a Markovian embedding by enlarging the dynamical space upon introduction 
of new state variables. The simplest two-state non-Markovian model with an 
exponentially decaying memory kernel can be embedded as a three state cyclic 
Markovian model, where the refractory period is exponentially distributed. Multi¬ 
state Markovian embedding also allows to treat a special Erlangian distribution of the 
refractory periods, which can be sharply peaked at a characteristic delay time. Moreover, 
models of bursting neurons and neurons with a power law distributed memory can, in 
principle, be considered in this generic Markovian embedding setup. The approach 
of Markovian embedding is especially suitable to treat the mean-field dynamics of the 
network, which presents a Markovian renewal process in the enlarged space of collective 
network variables. This is the simplest kind of network, where all the elements are 
virtually connected in all-to-all fashion. In this respect, the mean field dynamics of 
a network of non-Markovian renewal elements does not represent a renewal process in 
the reduced space of non-Markovian collective variables. Then, all the elements must 
be treated individually, keeping trace of their individual memory. The methodology of 
Markovian embedding allows to circumvent this problem for the mean-field dynamics. 

In the Wilson-Cowan model, two different sorts of interacting neurons are 
considered, excitatory and inhibitory. We focused on the simplest non-Markovian 
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generalization of this model, where the observed two-state non-Markovian dynamics of a 
single neuron is embedded as a three state cyclic Markovian process. The corresponding 
nonlinear mean-field dynamics is four dimensional. It has two dimensions more than in 
the original model. Moreover, it is stochastic and includes mesoscopic fluctuations due 
to a finite network size. For a sufficiently large system size, stochastic dynamics can be 
described within a Langevin equation approximation following the so-called chemical 
Langevin equation approach, with the noise terms vanishing in the deterministic limit 
of infinite size. We also exactly simulated the underlying dynamics as a continuous time 
Markovian random walk on a four-dimensional lattice using the well-known dynamical 
Monte Carlo (Gillespie) algorithm. The results of both stochastic approaches agree well 
with the deterministic dynamics within an oscillatory regime for a very large number 
of elements (several millions). Here, we showed that non-Markovian effects can be very 
essential. In particular, even deterministic dynamics with an exponentially decaying 
memory in the space of observable variables can be very different from the dynamics 
obtained within the Markovian approximation utilizing a delay-renormalized transfer 
function — the simplest approach to account for the delay effects. However, already 
the simplest approach allows to describe a dynamical phase transition from the silent 
network to coherent nonlinear oscillations of synchronized neurons upon a change of the 
delay period. This important feature is absent in the original Wilson-Cowan model. 

In more detail, we investigated the avalanche dynamics in a critically balanced 
network, where the processes of excitation and inhibition nearly compensate each 
other in the deterministic limit, where no avalanches are possible within the model 
considered. Mesoscopic noise fluctuations make, however, avalanches possible even in 
large networks with millions of neurons, where the deterministic description becomes 
completely inadequate, very differently from the oscillatory dynamics in such large 
networks. This result goes beyond the results in Ref. [15], where the avalanches cease 
to exist for already several tens of thousands elements. Even though a large delay 
should suppress avalanches by a transfer function renormalization if to think within 
the Markovian approximation, in reality the suppression is much weaker. Moreover, 
it turns out that the power law characterizing the distribution of avalanche sizes is 
very robust with respect to variation of both the delay period, and the system size, 
over several orders of magnitude, as well as the choice of the avalanche threshold. The 
latter fact proves that this is a real power law originated due to critical dynamics. It 
is characterized by a power-law exponent around a ~ —1.16 which is similar to the 
size exponent of the critical sandpile dynamics (though the both models are not really 
comparable). However, it is different from the critical exponent —1.62 found in Ref. [15], 
though for different network parameters. The distribution of the avalanches durations 
is, however, biexponential. We disproved that it presents a power law within our model, 
even though it can look like a power law over one time decade, as in experiments. In 
this respect, experiments mu seem to reveal a real power law with the critical size 
exponent a exp ~ —1.5 because its range extends with the growing system size. However, 
the experimental power law in the time duration does not show this important property. 
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As a matter of fact, it extends over merely one time decade being restricted by the period 
of 7 —oscillations. Any power law extending over one time or spatial decade can be fitted 
by a sum of just two exponentials, as we also show in this work for the time distribution. 
A further research is, therefore, required to clarify the nature of the apparent power law 
feature in the avalanche time distribution for the observed neuronal avalanches. 

Also very important is that the Langevin or diffusional approximation does change 
the critical exponents of the studied avalanche dynamics. There appears a power law in 
the time distribution, which is absent in the exact simulations, with the critical Langevin 
exponent 6 ^ ~ —1.06. Also the critical Langevin size exponent is different, ~ —1.03. 
This feature should be kept in mind while doing diffusional approximations in other 
models of critical dynamics. It may deliver incorrect results even for a large number of 
elements. 

We believe that the results of this work have methodological value and can be 
extended onto the dynamics of other networks with delay. They can serve also as a basis 
for further investigations of the role of non-Markovian memory effects in the dynamics 
of Wilson-Cowan type neuronal networks, including networks of bursting neurons, and 
networks with nontrivial topology, which we plan to investigate in future. 
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